################################################################################
# File name: figure_S3-S5.R
# Author: Dan M{\o}nster <danm@econ.au.dk>
# Date: 25 February 2024
# 
# Software versions used:
#     R Version 4.3.2 (2023-10-31)
#     dplyr Version 1.1.4
#     ggplot2 Version 3.5.0
#     patchwork Version 1.2.0
#
# Run this script with Rscript as: 'Rscript figure_S3-S5.R'
# Alternatively, run the script in an interactive R session, e.g. in RStudio.
################################################################################

# Load required packages using groundhog for reproducibility
#
# If you wish to run the script without the groundhog package, set the 
# use_groundhog variable below to FALSE.
#
# WARNING: This may affect the ability of this script to reproduce the results
# in the paper.
use_groundhog <- TRUE

if (use_groundhog) {
  options(repos = "https://cloud.r-project.org")
  library(groundhog)
  groundhog_day <- "2024-02-25"
  packages <- c("dplyr", "ggplot2", "patchwork")
  groundhog.library(packages, groundhog_day)
} else {
  warning("Not using groundhog. Script runs with available versions of packages.")
  library(dplyr)
  library(ggplot2)
  library(patchwork)
}


################################################################################
path <- "."
work_data_file <- paste0(path, "/../Input/Replicationdata.csv")
work_data <- read.csv(work_data_file)
out_path <- paste0(path, "/../Output/")
################################################################################


################################################################################
# Function used to create the plots in Figures S3, S4, and S5.
################################################################################
plot_ifsa_by_treatment <- function(data,
                                   figtitle  = "",
                                   ylabel = "Correct actions",
                                   figtag = "",
                                   jitter_w = 0.1,
                                   jitter_h = 0.2,
                                   bandwidth = 0.5,
                                   ybreaks = c(0, 2, 4, 6, 8, 10, 12),
                                   y_scale = 1) {
  # Summarise over recipes
  data <- data |> 
    group_by(treatment, pid) |> 
    summarise(IFSA = mean(IFSA),
              .groups = "drop")
  # Calculate mean and standard error in the two treatments
  data_mean_std <- data |> 
    group_by(treatment, pid) |> 
    summarise(IFSA = mean(IFSA),
              .groups = "drop") |> 
    ungroup() |> 
    group_by(treatment) |> 
    summarise(average = mean(IFSA),
              std_err = sd(IFSA) / sqrt(n()))
  # Create a data frame with number of data points
  annotation_data <- data |> 
    group_by(treatment) |> 
    summarise(N = n()) |> 
    mutate(text_label = paste0("N = ", N))
  # Set up a custom color palette
  treatment_colors <- c("Reminder" = "#4daf4a",
                        "ManyReminders" = "#377eb8",
                        "Control" = "grey")
  # Create the plot
  the_plot <- ggplot(data = data,
                     aes(x = treatment, y = IFSA / y_scale, fill = treatment)) +
    scale_color_manual(values = treatment_colors,
                       name = "treatment",
                       limits = unique(data$treatment),
                       aesthetics = c("colour", "fill")) +
    geom_violin(alpha = 0.7, bw = bandwidth) +
    geom_crossbar(data = data_mean_std |> 
                    mutate(average = average / y_scale,
                           std_err = std_err / y_scale),
                  aes(x = treatment,
                      y = average,
                      ymin = average - std_err,
                      ymax = average + std_err,
                      fill = treatment),
                  width = 0.9,
                  linewidth = 0.2,
                  fatten = 1) +
    geom_point(position = position_jitter(width = jitter_w, height = jitter_h),
               alpha = 0.4,
               colour = "black") +
    geom_text(data = annotation_data,
              aes(x = treatment, label = text_label),
              y = Inf) +
    # Avoid clipping geom_text label at y = Inf
    coord_cartesian(clip = "off") + 
    scale_y_continuous(breaks = ybreaks) +
    xlab("") +
    ylab(ylabel) +
    labs(tag = figtag, title = figtitle) +
    theme_classic() + 
    theme(legend.position = "none",
          plot.title = element_text(size = 10)) +
    expand_limits(y = 0)
  return(the_plot)
}


################################################################################
# Figure S3
################################################################################

# Common data used for the figure. Will be sub-set'ed for each panel below
data_S3 <- work_data |> 
  filter(treatment != 2,
         level == 8,
         tagrecipe != 0) |> 
  mutate(treatment = factor(treat_R,
                            levels = c(0, 1),
                            labels = c("Control", "Reminder")))

data_S3_A <- data_S3 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_score)

plot_S3_A <- plot_ifsa_by_treatment(
  data = data_S3_A,
  figtitle = "All actions\nFeedback study",
  figtag = "A",
  jitter_h = 0.3)
  
data_S3_B <- data_S3 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_remind)

plot_S3_B <- plot_ifsa_by_treatment(
  data = data_S3_B,
  figtitle = "Reminded actions\nFeedback study",
  figtag = "B",
  jitter_w = 0.2,
  jitter_h = 0.1,
  ybreaks = c(0, 1, 2, 3))

data_S3_C <- data_S3 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_not_remind)

plot_S3_C <- plot_ifsa_by_treatment(
  data = data_S3_C,
  figtitle = "Non-reminded actions\nFeedback study",
  figtag = "C",
  jitter_h = 0.3,
  ybreaks = seq(0,10,2))

data_S3_D <- data_S3 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_score)

plot_S3_D <- plot_ifsa_by_treatment(
  data = data_S3_D,
  figtitle = "All actions\nNo feedback study",
  figtag = "D",
  jitter_h = 0.3)

data_S3_E <- data_S3 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_remind)

plot_S3_E <- plot_ifsa_by_treatment(
  data = data_S3_E,
  figtitle = "Reminded actions\nNo feedback study",
  figtag = "E",
  jitter_w = 0.2,
  jitter_h = 0.1,
  ybreaks = c(0, 1, 2, 3))

data_S3_F <- data_S3 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_not_remind)

plot_S3_F <- plot_ifsa_by_treatment(
  data = data_S3_F,
  figtitle = "Non-reminded actions\nNo feedback study",
  figtag = "F",
  jitter_h = 0.3,
  ybreaks = seq(0,10,2))

# Arrange the plots into a figure using patchwork
figure_S3 <- (plot_S3_A + plot_S3_B + plot_S3_C) /
  (plot_S3_D + plot_S3_E + plot_S3_F)

# Save as PDF
S3_file <- paste0(out_path, "Figure_S3.pdf")
ggsave(S3_file, figure_S3, width = 4 * 3, height = 4 * 2)



################################################################################
# Figure S4
################################################################################

# Common data for Figure S4.
data_S4 <- work_data |> 
  filter(treatment != 0,
         level == 8,
         tagrecipe != 0) |> 
  mutate(treatment = factor(treat_MRR,
                            levels = c(0, 1),
                            labels = c("Reminder", "ManyReminders")))

data_S4_A <- data_S4 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_score)

plot_S4_A <- plot_ifsa_by_treatment(
  data = data_S4_A,
  figtitle = "All actions\nFeedback study",
  figtag = "A",
  jitter_h = 0.3,
  bandwidth = 0.7)
  
data_S4_B <- data_S4 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_remind)

plot_S4_B <- plot_ifsa_by_treatment(
  data = data_S4_B,
  figtitle = "Reminded in Reminder\nFeedback study",
  figtag = "B",
  jitter_w = 0.2,
  jitter_h = 0.1,
  ybreaks = c(0, 1, 2, 3),
  bandwidth = 0.7)

data_S4_C <- data_S4 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_many_remind)

plot_S4_C <- plot_ifsa_by_treatment(
  data = data_S4_C,
  figtitle = "Reminded in ManyReminders\nFeedback study",
  figtag = "C",
  jitter_w = 0.2,
  jitter_h = 0.1,
  ybreaks = seq(0,10,2),
  bandwidth = 0.7)

data_S4_D <- data_S4 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_never_remind)

plot_S4_D <- plot_ifsa_by_treatment(
  data = data_S4_D,
  figtitle = "Never reminded actions\nFeedback study",
  figtag = "D",
  jitter_w = 0.2,
  jitter_h = 0.2,
  bandwidth = 0.7)

data_S4_E <- data_S4 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_score)

plot_S4_E <- plot_ifsa_by_treatment(
  data = data_S4_E,
  figtitle = "All actions\nNo feedback study",
  figtag = "E",
  bandwidth = 0.7,
  jitter_h = 0.3)

data_S4_F <- data_S4 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_remind)

plot_S4_F <- plot_ifsa_by_treatment(
  data = data_S4_F,
  figtitle = "Reminded in Reminder\nNo feedback study",
  figtag = "F",
  jitter_w = 0.2,
  jitter_h = 0.1,
  ybreaks = c(0, 1, 2, 3),
  bandwidth = 0.7)

data_S4_G <- data_S4 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_many_remind)

plot_S4_G <- plot_ifsa_by_treatment(
  data = data_S4_G,
  figtitle = "Reminded in ManyReminders\nNo feedback study",
  figtag = "G",
  jitter_w = 0.2,
  jitter_h = 0.1,
  ybreaks = seq(0,10,2),
  bandwidth = 0.7)

data_S4_H <- data_S4 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_never_remind)

plot_S4_H <- plot_ifsa_by_treatment(
  data = data_S4_H,
  figtitle = "Never reminded actions\nNo feedback study",
  figtag = "H",
  jitter_w = 0.2,
  jitter_h = 0.2,
  bandwidth = 0.7)

# Arrange the plots into a figure using patchwork
figure_S4 <- (plot_S4_A + plot_S4_B + plot_S4_C + plot_S4_D +
    plot_layout(ncol = 4)) /
  (plot_S4_E + plot_S4_F + plot_S4_G + plot_S4_H +
     plot_layout(ncol = 4))

# Save as PDF
S4_file <- paste0(out_path, "Figure_S4.pdf")
ggsave(S4_file, figure_S4, width = 4 * 4, height = 4 * 2)


################################################################################
# Figure S5
################################################################################

# Common data for Figure S5
data_S5 <- work_data |> 
  filter(treatment != 2,
         level == 11,
         tagrecipe != 0) |> 
  mutate(treatment = factor(treat_R,
                            levels = c(0, 1),
                            labels = c("Control", "Reminder")))

data_S5_A <- data_S5 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_score)

plot_S5_A <- plot_ifsa_by_treatment(
  data = data_S5_A,
  figtitle = "All actions\nFeedback study",
  figtag = "A",
  jitter_h = 0.3)

data_S5_B <- data_S5 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_remind)

plot_S5_B <- plot_ifsa_by_treatment(
  data = data_S5_B,
  figtitle = "Reminded actions\nFeedback study",
  figtag = "B",
  jitter_h = 0.1,
  jitter_w = 0.2,
  ybreaks = c(0, 1, 2, 3))

data_S5_C <- data_S5 |> 
  filter(feedback == 1) |> 
  rename(IFSA = ifsa_not_remind)

plot_S5_C <- plot_ifsa_by_treatment(
  data = data_S5_C,
  figtitle = "Non-reminded actions\nFeedback study",
  figtag = "C",
  jitter_h = 0.3,
  ybreaks = seq(0,10,2))

data_S5_D <- data_S5 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_score)

plot_S5_D <- plot_ifsa_by_treatment(
  data = data_S5_D,
  figtitle = "All actions\nNo-feedback study",
  figtag = "D",
  jitter_h = 0.3)

data_S5_E <- data_S5 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_remind)

plot_S5_E <- plot_ifsa_by_treatment(
  data = data_S5_E,
  figtitle = "Reminded actions\nNo-feedback study",
  figtag = "E",
  jitter_h = 0.1,
  jitter_w = 0.2,
  ybreaks = c(0, 1, 2, 3))

data_S5_F <- data_S5 |> 
  filter(feedback == 0) |> 
  rename(IFSA = ifsa_not_remind)

plot_S5_F <- plot_ifsa_by_treatment(
  data = data_S5_F,
  figtitle = "Non-reminded actions\nNo-feedback study",
  figtag = "F",
  jitter_h = 0.3,
  ybreaks = seq(0,10,2))

# Arrange the plots into a figure using patchwork
figure_S5 <- (plot_S5_A + plot_S5_B + plot_S5_C) / (plot_S5_D + plot_S5_E + plot_S5_F)

# Save as PDF
S5_file <- paste0(out_path, "Figure_S5.pdf")
ggsave(S5_file, figure_S5, width = 4 * 3, height = 4 * 2)

# sessionInfo() # Run this after script has run to get all dependencies

